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Abstract 

We consider flame front propagation in channel geometries. The steady state 
solution in this problem is space dependent, and therefore the linear stability 
analysis is described by a partial integro-differential equation with a space 
dependent coefficient. Accordingly it involves complicated eigenfunctions. We 
show that the analysis can be performed to required detail using a finite order 
dynamical system in terms of the dynamics of singularities in the complex 
plane, yielding detailed understanding of the physics of the eigenfunctions 
and eigenvalues. 

1. INTRODUCTION 

In this paper we discuss the stability of steady flame fronts in channel geometry. Tradi- 
tionally [1-3] one studies stability by considering the linear operator which is obtained by 
linearizing the equations of motion around the steady solution. The eigenfunctions obtained 
are delocalized and in certain cases are not easy to interpret. In the case of flame fronts the 
steady state solution is space dependent and therefore the eigenfunctions are very different 
from simple Fourier modes. We show in this paper that a good understanding of the nature 
of the eigenspectrum and eigenmodes can be obtained by doing almost the opposite of tradi- 
tional stability analysis, i.e., studying the localized dynamics of singularities in the complex 
plane. By reducing the stability analysis to a study of a flnite dimensional dynamical system 
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one can gain considerable intuitive understanding of the nature of the stabihty problem. 

The analysis is based on the understanding that for a given channel width L the steady 
state solution for the flame front is given in terms of N{L) poles that are organized on a line 
parallel to the imaginary axis [4]. Stability of this solution can then be considered in two 
steps. In the first step we examine the response of this set of N{L) poles to perturbations 
in their positions. This procedure yields an important part of the stability spectrum. In the 
second step we examine general perturbations, which can also be described by the addition 
of extra poles to the system of N{L) poles. The response to these perturbations gives us 
the rest of the stability spectrum; the combinations of these two steps rationalizes all the 
qualitative features found by traditional stability analysis. 

In Sec. 2 we present a brief review of the stationary solutions of front propagation in 
channel geometries. In Sec. 3 we present the results of traditional linear stability analysis, 
and show the eigenvalues and eigenfunctions that we want to interpret by using the pole 
decomposition. Sec. 4 presents the analysis in terms of complex singularities, in two steps 
as discussed above. A summary and discussion is presented in Sec. 5. 



2. FLAME PROPAGATION IN CHANNEL GEOMETRY 

We consider a channel of transverse width 27rL, and of infinite extent in the longitudinal 
direction. A graph of a fiame front propagating in this channel is defined as h{9, t) where 9 is 
a rescaled transversal coordinate < ^ < 27r. The equation of motion is given conveniently 
in terms of u{e, t) = dh(e, t)/de [5-8]: 

du{e,t) _ u{e,t)du{e,t) v d'^u{e,t) 

+\i{u{e,t)], (1) 

Here is a viscosity- like parameter and the functional I[u{6,t)] is conveniently defined in 
terms of the spatial Fourier transform 

u{e,t) = r e'^'^u{k,t)dk , (2) 

J —oo 
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I[u{k,t)] = \k\u{k,t) . 



(3) 



It is very useful [4,9-11] to discuss the solutions of these equations of motion in terms of 
expansions in N poles whose position Zj{t) = Xj{t) + iyj{t) in the complex plane is time 
dependent: 



N 

u{9, t) = cot 

N 



9-zj{t) 



+ c.c. 



2sin[^ - Xj{t)] 



(4) 



cosh[yj{t)] - cos[9 - Xj{t)] 

Substituting (4) in (1) we derive the following ordinary differential equations for the 
positions of the poles: 



dx- ^ 

k=i,k=ij 



cosh(yj - yk) 



(5) 



cos{xj — Xk) + cosh.{yj + yk) — cos{xj — Xk) 

sinh(|/j - yk) 

' ^cosh(|/j- - yk) - cos{xj - Xk) 

smh{yj + yk) \ , . 

r + u cothfy,) - L. 



dt 



(6) 



cosh(|/j + yk) - cos{xj - Xk) 

In particular we can find the steady state solution Us{0) by demanding Xj — yj — and 
stabihty. The solution is 



TV 



u 



2 sin[^ — Xs 



(7) 



cosh[yj] - cos[6' - Xs] 

where Xg is the real (common) position of the stationary poles and yj their stationary 
imaginary position. We need to determine the actual positions yj. This is done numerically 
by running the equations of motion for the poles starting from N poles in initial positions and 
waiting for relaxation. A complete analysis of this steady-state solution was first presented 
in Ref. [4] and the main results are summarized as follows: 

1. There is only one stable stationary solution which is geometrically represented by 
a giant cusp (or equivalently one finger) and analytically by N{L) poles which are 



aligned on one line parallel to the imaginary axis. The existence of this solution is 
made clearer with the following remarks. 

2. There exists an attraction between the poles along the real line. This is obvious from 
Eq.(5) in which the sign of dxj/dt is always determined by sin(a;j — Xk). The resulting 
dynamics merges all the x positions of poles whose y-position remains finite. 

3. The y positions are distinct, and the poles are aligned above each others in positions 
Uj-i < Uj < Uj+i with the maximal being yN{L)- This can be understood from Eq.(6) 
in which the interaction is seen to be repulsive at short ranges, but changes sign at 
longer ranges. 

4. If one adds an additional pole to such a solution, this pole (or another) will be pushed 
to infinity along the imaginary axis. If the system has less than N{L) poles it is 
unstable to the addition of poles, and any noise will drive the system towards this 
unique state. The number N{L) is 

NiL) =\U- + ^)]^ (8) 



2 \u 

where ... is the integer part. To see this consider a system with N poles and such 
that all the values of yj satisfy the condition < yj < ymax- Add now one additional 
pole whose coordinates are Za = {xa, ya) with ya ^ y-max- From the equation of motion 
for ya-i (6) we see that the terms in the sum are all of the order of unity as is also the 
cot(ya) term. Thus the equation of motion of ya is approximately 

dya 2Ar + 1 1 



dt L2 L 



(9) 



The fate of this pole depends on the number of other poles. If is too large the pole 
will run to infinity, whereas if N is small the pole will be attracted towards the real 
axis. The condition for moving away to infinity is that N > N{L) where N(L) is given 
by (8). On the other hand the y coordinate of the poles cannot hit zero. Zero is a 
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repulsive line, and poles are pushed away from zero with infinite velocity. To see this 
consider a pole whose Uj approaches zero. For any finite L the term coth(|/j) grows 
unboundedly whereas all the other terms in Eq.(6) remain bounded. 

5. The height of the cusp is proportional to L. The distribution of positions of the poles 
along the line of constant x was worked out in [4] . 

We will refer to the solution with all these properties as the Thual-Frisch-Henon (TFH)-cusp 
solution. 



3. LINEAR STABILITY ANALYSIS IN CHANNEL GEOMETRY 

The standard technique to study the linear stability of the steady solution is to perturb 
it by a small perturbation (f){9,t): u{9,t) — Us{0) + (f){9,t) . Linearizing the dynamics for 
small (j) results in the following equation of motion 

d(i>{e,t) _ 1 r 



dt L2 



de[us{e)(i){e,t)] 



+ ud'o(i>{9,t)] + ^i{(l>{e,t)) . (10) 

were the linear operator contains Us{9) as a coefficient. Accordingly simple Fourier modes do 
not diagonalize it. Nevertheless, we proceed to decompose (f){x) in Fourier modes according 

to , 

oo 

mt)= E Mt)e''' (11) 

k=—oo 

oo JV 

Us{9) = -2ui E E5^^^('^)e"'*''^'e'*'^ (12) 

fc= — CO j = l 

The last equation follows from (7) by expanding in a series of sin k9. In these sums the 
discrete k values run over all the integers. Substituting in Eq.(lO) we get: 

dkjt)) 

dt 

where akn are entires of an infinite matrix: 



^akn(i>n{t) , (13) 



akk^ -j^- j^k^ , (14) 
akn = ^szgn{k-n){2iyJ2e-^''-''^''^) k^n. (15) 

To solve for the eigenvalues of this matrix we need to truncate it at some cutoff /c- vector k*. 
The scale k* can be chosen on the basis of Eq.(14) from which wc sec that the largest value 
of k for which Ukk > is a scale that we denote as k^ax, which is the integer part oi L/v. 
We must choose k* > k^ax and test the choice by the convergence of the eigenvalues. The 
chosen value of k* in our numerics was 4:kmax- One should notice that this cutoff limits the 
number of eigenvalues, which should be infinite. However the lower eigenvalues will be well 
represented. The results for the low order eigenvalues of the matrix akn that were obtained 
from the converged numerical calculation are presented in Fig.l 

The eigenvalues are multiplied by L'^/i/ and are plotted as a function of L. We order the 

eigenvalues in decreasing order and denote them as Aq > Ai > A2 The figure offers a 

number of qualitative observations: 

1. There exists an obvious Goldstone or translational mode u'g{6) with eigenvalue Aq = 0. 
This eigenmode stems from the Galilean invariance of the equation of motion. 

2. The eigenvalues oscillate periodically between values that are L-independent in this 
presentation (in which we multiply by L^). In other words, up to the oscillatory 
behaviour the eigenvalues depend on L like L"^. 

3. The eigenvalues Ai and A2 hit zero periodically. The functional dependence in this 
presentation appears almost piece-wise linear. 

4. The higher eigenvalues are more negative. They exhibit similar qualitative behaviour, 
but without reaching zero. We note that the solution becomes marginally stable for 
every value of L for which the eigenvalues Ai and A2 hit zero. The dependence of 
the spectrum indicates that the solution becomes more and more sensitive to noise as 
L increases [12]. 
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In addition to the eigenvalues, the truncated matrix also yields eigenvectors that we denote 
as A^^\ Each such vector has k* entries, and we can compute the eigenfunctions f''^\6) of 
the linear operator (10), using (11), as 

Eq.(lO) does not mix even with odd solutions in 6, as can be checked by inspection. Con- 
sequently the available solutions have even or odd parity, expandable in either cos or sin 
functions. The first two nontrivial eigenfunctions f^^\0) and f^'^\9) are shown in Figs. 2, 3. 
It is evident that the function in Fig.2 is odd around zero whereas in Fig. 3 it is even. Sim- 
ilarly we can numerically generate any other eigenf unction of the linear operator, but we 
understand neither the physical significance of these eigenfunction nor the L dependence of 
their associated eigenvalues shown in Fig.l In the next section we will demonstrate how the 
dynamical system approach in terms of singularities in the complex plane provides us with 
considerable intuition about these issues. 

4. LINEAR STABILITY IN TERMS OF COMPLEX SINGULARITIES 

Since the partial differential equation is continuous there is an infinite number of modes. 
To understand this in terms of pole dynamics we consider the problem in two steps: First, we 
consider the 2N[L) modes associated with the dynamics of the N{L) poles of the giant cusp. 
In the second step we explain that all the additional modes result from the introduction of 
additional poles, including the reaction of the N{L) poles of the giant cusp to the new poles. 
After these two steps wc will be able to identify all the linear modes that were found by 
diagonalizing the stability matrix in the previous section. 

A. The modes associated with the giant cusp 

In the steady solution all the poles occupy stable equibirilium positions. The forces 
operating on any given pole cancel exactly, and we can write matrix equations for small 
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perturbations in the pole positions 5yi and Sxi. 

Following [4] we rewrite the equations of motion (6) using the Lyapunov function U : 



where i — 1, N and 



U = — In sinh yi + 2 ^(In sinh 



i<k 



+ lnsinh^^)]-^y, (18) 
The linearized equations of motion for 5yi are: 

The matrix d'^U/dyidyk is real and symmetric of rank N. We thus expect to find N real 
eigenvalues and N orthogonal eigenvectors. 

For the deviations Sxi in the x positions we find the following linearized equations of 
motion 

ly ^ 1 

L5xj^—-5xj Yl 



L ^j^^^^. cosh(y,- - Vk) - 1; 

) 



cosh(|/j +yk)-l' 
u ^ 1 1 

fe=lX7^i cosh(yj - yk) - 1) cosh(yj- + y^) - 1 



In shorthand: 



= . (21) 

The matrix V is also real and symmetric. Thus V and d'^U / dyidyk together supply 2N{L) 
real eigenvalues and 2N{L) orthogonal eigenvectors. The explicit form of the matrices V 
and d'^U / dyidyk is as follows: For i k: 

d'^U u. 1/2 1/2 ^ 



dyidyk L^smh\^) sinh2(^ 



(22) 



L cosh{yi -yk)-l) cosh(|/j + yk) - 1 



and ioT i — k one gets: 



T 



dyf L sinh'(^) 2 J 

"^sinh^(l/i)^ ^^^^ 

^ z/ 1 1 

" §5^~L(cosh(y,-yfe)-l) ^ cosh(y, + yfc)-l^^ ^^^^ 

Using the known steady state solutions yi at any given L we can diagonalize the N{L) x N{L) 

matrices numerically. In Fig.4 we present the eigenvalues of the lowest order modes obtained 

from this procedure. The least negative eigenvalues touch zero periodically. This eigenvalue 

can be fully identified with the motion of the highest pole ywiL) in the giant cusp. At isolated 

values of L the position of this pole tends to infinity, and then the row and the column in our 

matrices that contain yN{L) vanish identically, leading to a zero eigenvalue. The rest of the 

upper eigenvalues match perfectly with half of the observed eigenvalues in Fig.l. In other 

words, the eigenvalues observed here agree perfectly with the ones plotted in this Fig.l until 

the discontinuous increase from their minimal points. The "second half" of the oscillation 

in the eigenvalues as a function of L is not contained in this spectrum of the N{L) poles of 

the giant cusp. To understand the rest of the spectrum we need to consider perturbation of 

the giant cusp by additional poles. The eigenfunctions can be found using the knowledge of 

the eigenvectors of these matrices. Let us denote the eigenvectors of d'^U/dyidyk and V as 

a^^^ and b^^^ respectively. The perturbed solution is explicitly given as (taken for Xs — 0): 

,- sm(9 — 6xi) 

^ - g coshfa + 4.)-oos(.-.x.) 

where 5m is 

N oo 

5u^ - 4i/ ^ ^ Syike'''^' sin k9 

i=l k=l 

N oo 

-4:uJ2J2 Sxike-'^y' cos kO (27) 

1=1 k=l 



So knowing the eigenvectors a*^^) and h^^^ we can estimate the eigenvectors f^^\9) of (16): 

N oo 

/iV) = -4^^ E E "F^e-'^^^ sin k9, j = 1, N (28) 

i=l A;=l 

or 

N oo 

(e) = -41/ E E ^'?^A;e-'=^^ cos kO , j^l,...,N (29) 

i=l fc=l 

where we display separately the sin expansion and the cos expansion. For the case j — 1, 
the eigenvalue is zero, and a uniform translation of the poles in any amount Sxi results in a 
Goldstone mode. This is characterized by an eigenvector bf^^ = 1 for all i. The eigenvectors 
fW (Fig.5,6)computed this way are identical to numerical precision with those shown in 
Figs. 2, 3, and observe the agreement. 



B. Modes related to additional poles 

In this subsection we identify the rest of the modes that were not found in the previous 
subsection. To this aim we study the response of the TFH solution to the introduction 
of additional poles. We choose to add M new poles all positioned at the same imaginary 
coordinate yp <^ Umax, distributed at equidistant real positions {xj — xq + (27r/M)_7}^^. 
For xo = we use (4) and the Fourier expansion to obtain a perturbation of the form 

5u{e, t) ~ 4i/Me"^^^(*) sin MO (30) 

For xo — — 7r/2M we get 

6u{9, t) ~ 4i/Me-^^^(*) cos M9 (31) 
in both cases the equations for the dynamics of follow from Eqs.(5)-(6): 

f^2ii«(M), (32) 

where a{M) is given as: 
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a(M) = [l(^ + l)]-l(^-M) (33) 

Since (32) is linear, we can solve it and substitute in Eqs.(30)-(31). Seeking a form 5u{d, t) ~ 
exp{—\{M)t) we find that the eigenvalue A(M) is 

MM) = 2M-^a(M) (34) 

These eigenvalues are plotted in Fig. 7 At this point we consider the dynamics of the poles 
in the giant cusp under the influence of the additional M poles. From Eqs.(19), (21), (5), 
(6) we obtain, after some obvious algebra, 

LSy^ = E ^^^y^ - 4yMe-^^^W sinh(My,) (35) 

or 

L5xi = ^ij^^j - 4yMe-^^^(*) cosh(Myi) (36) 

It is convenient now to transform from the basis Syi to the natural basis Wi which is obtained 
using the linear transformation w — A~^Sy. Here the matrix A has columns which are 

the eigenvectors of d'^U/dyidyj which were computed before. Since the matrix was real 
symmetric, the matrix A is orthogonal, and A~^ — A^. Define C — A-^Me'^'^^^^^ and write 

Wi = -XiWi - Ce-^(^)*ei , (37) 
where — Aj are the eigenvalues associated with the columns of A, and 

= E^iiSinhMy,- . (38) 

i 

We are looking now for a solution that decays exponentially at the rate A(M): 

Wi{t) = Wi(0)e-^W* (39) 
Substituting the desired solution in (37) we find a condition on the initial value of wf 



Transforming back to 5yi we get 

C 

%(0) = E AfeWfc(O) = - E AfcT rTTTT E ^'fe sinh My, 

= -CEsmhMy,El^ (41) 

I k 

We can get the eigenfunctions of the hnear operator, as before, using Eqs.(27), (30), (31), 
(41). We get 

. X ^^^^ °° A A 

j=l k=l I m ^m- A(^MJ 

X ke-'^y' sin A;^ + L^C sin MO (42) 

An identical calculation to the one started with Eq. (37) can be followed for the deviations 
Sxi. The final result reads 

NiL) oo 4 4 

f^^\0)^ACu E E(EcoshM,/^E r ) 

i=l k=l I m ^m — ^[M) 

xke-^y^ COS ke + L^C cos MO , (43) 

where A is the matrix whose columns are the eigenvectors of V , and — Aj its eigenvalues. 

We are now in position to explain the entire linear spectrum using the knowledge that 
we have gained. The spectrum consists of two separate types of contributions. The first 
type has 2N modes that belong to the dynamics of the unperturbed iV(L) poles in the giant 
cusp. The second part, which is most of the spectrum, is built from modes of the second 
type since M can go to infinity. This structure is seen in the Fig. 4 and Fig. 7. 

We can argue that the set of eigenfunctions obtained above is complete and exhaustive. 
To do this we show that any arbitrary periodic function of can be expanded in terms 
of these eigenfunctions. Start with the standard Fourier series in terms of sin and cos 
functions. At this point solve for sinkO and cos kO from Eqs. (42-43). Substitute the results 
in the Fourier sums. We now have an expansion in terms of the eigenmodes /(^) and in 
terms of the triple sums. The triple sums however can be expanded, using Eqs. (28-29), 
in terms of the eigenfunctions /*^^\ We can thus decompose any function in terms of the 
eigenfunctions /(^) and /^^-*. 
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5. CONCLUSIONS 



We discussed the stability of flame fronts in channel geometry using the representation of 
the solutions in terms of singularities in the complex plane. In this language the stationary 
solution, which is a giant cusp in configuration space, is represented by N{L) poles which are 
organized on a line parallel to the imaginary axis. We showed that the stability problem can 
be understood in terms of two types of perturbations. The first type is a perturbation in the 
positions of the poles that make up the giant cusp. The longitudinal motions of the poles 
give rise to odd modes, whereas the transverse motions to even modes. The eigenvalues 
associated with these modes are eigenvalues of a finite, real and symmetric matrices, cf. 
Eqs.(22), (23), (24), (25). The second type of perturbations is obtained by adding poles to 
the set of N(L) poles representing the giant cusp. The reaction of the latter poles is again 
separated into odd and even functions as can be seen from Eqs.(30), (31). Together the 
two types of perturbations rationalize and explain all the features of the eigenvalues and 
eigenfunctions obtained from the standard linear stability analysis. 
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FIG. 1. A plot of the first five eigenvalues obtained by diagonalizing the matrix obtained by 
traditional stability analysis, against the system size. The eigenvalues are normalized by I? jv. 
The largest eigenvalue is zero, which is a Goldstone mode. All the other eigenvalues are negative 
except for the second and third that touch zero periodically. The second and fourth eigenvalues 
are represented by a solid line and the third and fifth eigenvalues are represented by a dot-dashed 
line. 
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FIG. 2. The first odd eigenfunction obtained from traditional stability analysis. 
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FIG. 3. The first even eigenfunction obtained from traditional stability analysis. 
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FIG. 4. The eigenvalues associated with perturbing the positions of the poles that consist the 
giant cusp. The largest eigenvalue is zero. The second, third, fourth and fifth eigenvalues are 
represented by a solid line, dot-dashed line, dotted line and dashed line respectively. 
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FIG. 5. The first odd eigenfunction associated with perturbing the positions of the poles in the 
giant cusp. 
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FIG. 6. The first even eigenfunction associated with perturbing the positions of the poles in 
the giant cusp. 
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FIG. 7. Spectrum of eigenvalues associated with the reaction of the poles in the giant cusp to 
the addition of new poles. 
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